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1. One~dlinenslonal Splines 


The purpose of this section Is to give an Introduction to some outstanding 
properties of one -dimensional splines. We are not interested In giving a very 
detailed discussiwi of this vast subject but only in pointing out some results 
which serve to Illustrate some properties of multivariate splines. 

1.1. The Cubic Spline 

Consider knots ^2 ^ ^3 *^ * * * "^ ^'n ® real interval [0,T] and 

real data y^^, •••> yjj* oldest problems considered by mathema- 

ticians and engineers has been that of finding a "smooth" curve joining the 

2 

points ^ IR . Or, in mathematical terms, to find a function 

g: [0,T] -*-lR such that g takes the value y^ at t^ i = l,...,n. 

g(t^) = l,...,n (1.1.1) 

Obviously, this problem does not have a unique solution and the choice 
of a solution will depend on many factors, some of which are: 

— Smoothness properties 
— Convergence properties 
— Cost of computation. 

Even if the first two properties have a great importance, the third one is 

usually the one determining the interpolant to be used. Among the solutions 

# 

used in the past, the most popular one. was the polynomial interpolation. As 
it is known, we can find a unique polynomial p^ of degree n-1 satisfying 
(1.1.1). Moreover, an explicit expression for p^ can be given in terms of 
Lagrange polynomials 
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where 


n 

n 

i«i 




£^(x) 


n (x-x ) 

iii L 

w' (x^) 


( 1 . 1 . 2 ) 


a .1.3) 


n 

w(x) = n (x-x.) . (1.1.4) 

i=l 

Thus, the solution can be easily written as a linear combination 

of well-known polynomials. This is a very interesting property, even if 
(1.1.2) is not used for practical calculations because it is cheaper to use 
the Newton's formula (cf. [ 6 ]) . 

Another interesting property of the interpolation polynomial is that it 
can be differentiated infinitely many times and each derivative is still a 
polynomial and can be easily evaluated using a Homer's scheme. This last 
assertion is only partly true because the evaluation of a polynomial of 
degree n-1 requires n-1 multiplications and for large n this could become 
very expensive. 

However, the main trouble arising with pol 3 momall interpolation is mainly 
due to the fact that they arc too smooth and hence too inflexible. This leads 
to oscillations when the number of points is high (more than 20) . Indeed, 
the interpolation polynomials have very bad approximation properties except 
for a small class of analytic functions (cf. [ 6 ]). 

Schoemberg [ 25 ] found the way to avoid the problems coming from the 

Inflexibility of polynomials and still keep their nice properties; ht intro- 


duced the use of piecewise polynomials matched together by continuity condi- 
tions. He called the new class of functions: splines. 

For a general description of polynomial splines the reader is referred 
to [ 1 ], [ 9 ], [ 26 ]. In this section we will only describe one 

of the most popular splines: the cubic spline. 

Let be the set of all functions satisfying the following 

properties: 

CSl: a is a cubic polynomial in i = l,...,n-l 

CS2: o e C^(t, ,t ) 

1 n 

It is well known that is a linear space of dimension n+2. So» if 

we want to use as interpolant a function a e we mu^t impose two additional 
conditions in order to determine the n+2 free parameters. Among the most 
popular conditions we have: 


I. Natural 

a”(t^) 


11. Periodic 

a’(t^) 

- 


a'*(t^) 

■ 

III. Hermite 

a’(t^) 



a-'(tn) 

- • 


Optimal convergence rates are obtained using conditions type II or III. 
In these lectures we will deal only with natural conditions . For a more detailed 
description of convergence rates the reader is referred to [ 9 ], [ 26 ] 


and the references therein. 


Let be the linear space of natural cubic splines, that Is, the 

space of elements of satisfying condition I. It Is an n-dlmenslonal 

linear space and for any y^^, ® there exist one and only one 

element In such that 

a(ti) = i “ 1 n • 

This new Interpolant keeps some of the most Important properties of 
polynomials : 

1) The evaluation of o at a given point t only requires multiplica- 
tions and additions 

2) The integral or derivative of a is still a piecewise polynomial. 

The advantage of this kind of interpolant is that we get "good" conver- 
gence rates even for functions being only continuous. (For details, see 

[ 9 ], [26 ]). 

Of course we do not have an explicit formula like (1.1.2) and we will have 
to solve a linear system in order to find the n free parameters determining 
a. Many methods are available, but we would like to point out one that is 
now able to give a generalization of splines to two and higher dimen- 
sions. It is the B-spline method. (Here B stands for basis) (cf. [ 5 ], 

I 9 ], [26 ]). 

The idea is to construct a basis for with local support, that is, 

such that each element in the basis has the smallest possible support. Then 
the evaluation of o at a point t needs the calculation of a linear combina- 
tion with only a small number of nonzero terms. 

More precisely, let x _2 - '** “ *** ^ *n+3 extended 

partition 
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^-2 ■ ^-1 " *0 " h 


Xi . 


1 ■ 1, . . . ,n 


X., ■ x._ ■ Xj^, “ t 

n+1 n+2 n+3 n 


Now define the nf2 B-spllnes by the formula 


*l<t) • 


i-2’*i-l*^l*^i+l‘^i+2 


(t-x)^ , 


(1.1.5) 


t, £ t < t , i ■ 0, ...,n+l , 

1 n ’ » » * 


where 5 f is the divided difference of f with respect to 

Ui- « • • • * uc 

!■ P 

3 

a . And the truncated power function (t-x)_^ is given by 


(t-x) 


i (t-x) 


t < X 

t 2: X 


( 1 . 1 . 6 ) 


For t “ t we have 
n 




lim B (t) 

t-*-t“ 

n 


(1.1.7) 


The set {Bq, . . . is a basis for with the following interesting 

properties; 

T heorem 1.1.1 . The B-splines B^, defined by (1.1.5) - (1.1.7) 

satisfy 

( 1 . 1 . 8 ) 


Bj(t) - 0 




t « ^^i_2»^i+2^ ' 


B^(t) > 


0 


(1.1.9) 


r 
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Moreover, for any t « have 


n+1 

J B (t) - 1 . (1.1.10) 

1-0 ^ 


Proof. See [ 26 ]. 


These remarkable properties imply among others that the linear system to 
be solved has a band matrix. More precisely, the matrix Is trldlagonal In 
the cubic case. Another Interesting fact extensively used for computational 
purposes Is the recurrence relation which allows a very simple algorithm for 
the calculation of B^(t) . For details see [ 26 ] . 

As we shall see In other sections, the first extension of spline techniques 
to higher dimensions has been given using tensor products of one-dlmenslonal 
splines. The small support basis for that case is easily found using tensor 
products of B-spllnes. This kind of technique can be used in Interpolation, 
smoothing or curve fitting, as we shall see later. 

More recently there has been a lot of Interest in the development of multi- 
dlmens:'.or.al B-splines. We are not going to discuss this approach during these 
lectures and the interested reader Is referred to [ 7 ], [ 8 ] and 

the references therein. 

Now we turn our attention to a different property of polynomial splines 
which allows a very natural extension to higher dimensions. We begin with a 
lemma. 

2 

Lemma 1.1.2 . Let g « C [tj^,t^], then for any b have: 


g”(t)s"(t) dt 


n 

I (s 

1-1 


’"(t^) - s"'(t^)) g(t^) 


4 


( 1 . 1 . 11 ) 
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where 

8”'(t") E 0 

s”'(t") = 0 . 

n 

Proof . Integrating by partt^, using the fact that s'" Is piecewise constant 

and rearranging terms we get the desired result. 

2 

Now consider any function y in C interpolating the data 

^1’ ’*’* ^n ^1’ ^n* ^ unique element in inter- 

polating the data. Then from (1.1.11) we have 

t^ 

/ (y"(t) - o"(t)) a"(t) dt - 0 . 

^1 

Or, if we extend a to [0,T] by 

o(t^) + (t - t^)a' (tj^) t < t, 

o(t) * a(t) t, < t < t 

1 n 

a(t ) + (t - t )a’(t ) t > t 

'■ n n n n 

we have 

T 

/ (u"(t) - o"(t)) a"(t) dt ■ 0 . (1.1.12) 

0 

More precisely, we have the following: 

2 

Theorem 1.1.3 . For any y « H [0,T] interpolating y^^ y^ at t^ , ..., t^, 

we have 

(y-o,a) ■ 0 

where o is the natural cubic spline interpolating the data and 
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T 

(y.v) - / y"(t) v”(t) dt . (1.1.13) 

0 

Proof . See [ 19 ], I 26 ]. 

From this theorem, usually called the First Integral Relation, we Immedia- 
tely get the following property that Is usually tak&i as an equivalent defini- 
tion of the natural cubic spline: 

Theorem 1.1.4 . The natural cubic spline o interpolating y^^, ..., y^ 
at t^, ..., y^ Is the unique solution to the following minimization problem. 

Minimize (y,y) (1.1.14) 

y e I 

y 

where 

ly - (y « H^[0,T] 1 y(t^) - y^^, i “ l,...,n> . 

Proof . From the First Integral Relation we get 


(y,y) 


(y-o+o, y-o-Kj) 


(y-o, 

(y-o. 


u 

y-o) + 2(oyy-o) 
y-o) + (o,o) > 


+ (o,o) 

(o,o) , 


y e I 

y 


This proves that o minimizes (*,*) over 1^. The uniqueness comes 

from the fact that the kernel of the semi-norm (*,*) is P^^ the set of 

2 

polynomials of degree one whose intersection with 1^ ■ (y « H [0,T] | 

Is the element 0. This concludes the proof. 
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This vsrlscional definition of the natural cubic splines has lead to 
several generalizations in one and many dimensions. We are going to study 
some of these during these lectures. For additional details see [ 3 ], 

[ 7 ], [ 9 ], [ 13 ], [ 19 [ 26 ] end the references therein. 


1.2. Smoothing Splines 

Consider now the case tdiere the data are noisy. This is the case in most 
practical problems. Assume the model 

y^^ - f(t^) + i ■ l,«..tn (1.2.1) 

where f is the "smooth" unknown function to be approximated and the errors 
i ■■ l,...,n are assumed to be realizations of independent identically 
distributed normal random variables with 


E[e^l - 0 

E[e^] - i - l,...,n 

E[e^£j] - 0 1 J . 


To solve this problem Schoemberg introduced the smoothing spline. This 
spline is the unique solution to the following minimization problem. 


Minimize 

P«H^10,T] 


■X/\u"(t))^dt + ^ [ (y(t,)-y,)4 
0 "l-l ^ ^ j 


( 1 . 2 . 2 ) 


Where 1 > 0 is the smoothing parameter controlling the tradeoff between 

the smoothness of the solution measured by (y,y) and the approximation to 

1 " 2 
the data measured by — ^ (y(t. ) - y^)^. 

" i-1 ^ 
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The solution to (1.2.2) Is sn element of that is, a natural cubic 


spline. Moreover, let ^ be the nxn symmetric semi-positive definite matrix 


such that 


2^ Z 


min (y,M) 
U«H^[0,T] 


W(ti)-2i 


i"i, . . . ,n 


This matrix is well defined because the transformation associating the 


Interpolating cubic spline to the data is linear and (*,*) is 


quadratic. 


With this definition (1.2.2) can be restated as: 


( f in 

Minimize ^Minimize lX(y,vi) + — (u(t ) -y.)‘ 

z.iP 1 




1* 1 1 . . . , n 


T 1 

Minimize • \ z Qz + — Y (y - z.) 

n. , 1 1 

z.lf 1 


And the solution z to this problem is then given by 


z ■ (I + nXf2) ■‘■y 


If we call 0 . to the smoothing spline then z is the vector of the 

XI ^ A 


values cf a . at the knots 
n,X 


^ ■ ‘’n.x'V 


Vi/ 


1 - 1 ,..., 


(1.2.6) 



« V i 
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As we shall see later, (1.2.5) Is not of practical use but it allows us 
to give a nice Interpretation of the effect of this procedure on the data. 

il being symnetric and positive semi-def Inlte, Its eigenvalues are positive 
real numbers. Moreover, we have (cf. [11 ]) for 1 ■ l,...,n; the 

eigenvalues of T. 

0 “ - P 2 < P 3 < ••• < Pjj . (1.2.7) 

Also, for n sufficiently large and 1 > 10 (for practical purposes), we 
have (cf. [ 33 ]), 

i** 

p^ a — (1.2.8) 

1 n 

Let Q be the orthogonal matrix diagonalizing then 

= Q 2 Q* (1.2.9) 

where Z = diag (p, , p_, ..., p ). And 
1 z n 

z = Q(I + nXZ) ^Q*y 

or 

z = (I + nXE)~^y (1.2.10) 

with 

z = Q*z 

y = Q*y . 

If we write (1.2.10) In terms of each component and use (1.2.8) we obtain 

"i iT^yi • 


(1.2.11) 


Thus we are filtering the data using a lew-pass Butterworth type filter 

-1 

with cutting frequency (Xa)’'^ to be chosen. 

The choice of X from the data Is a very complicated problem that has 
received a lot of attention during the last years. For details see [ 4 ], 

[ 17 ], [ 30 ], [ 31 ], I 37 ], [ 38 ] . We are just going to describe 

two of them. The first one was proposed by Relnsch [ 24 ] and is closely 

related to the methods for choosing the regularization parameter in Tikhonov 
Regularization [ 38 ] . 

The idea is to solve the following equation for 


or, equivalently 




= a 


n n^X^pf 


1 Y ^i -v2 

n (1+nXp^) ^ ^i 


( 1 . 2 . 12 ) 


an equation that can be solved using bisection or another nesting method. 

As it can be seen, the method is indeed very simple and the idea is 
attractive. Nevertheless thert is a lot of numerical and theoretical evidence 
showing that in most cases the solution obtained by this method gives a spline 
which tends to be too smooth, i.e., it tends to eliminate variations of the 
true function (cf. [44 ]) . 

The second method is called Generalized Cross-validation and was introduced 
in this form by Craven and Wabba [ 4 ]. The idea of the method is more diffi- 
cult to introduce and the interested reader should go to their paper [ 4 ] 
for a detailed description. See also [ 17 ], [ 30 ] , [ 31 ] . 


The GCV method reduces to the minimization over X of the following loss 
function, usually referred to as the GCV function: 


where 


V(X) 


i*l ’ 


fl - ^tr.(A(X))l^ 
n I 


A(X) = (1 + nXQ) 


(1.2.13) 


Special methods have been developed for the efficient con^jutation of the 
minimi 2 er of V. For details see [ 30 ] , [ 31 ] , [ 43 ] . 


2. Thin Plate Splines 

In this section we introduce one of the most popular forms of bldimenslonal 
splines and give theoretical and computational details. We have chosen to start 
with these kind of splines because we think it is a natural generalization of 
the popular cubic spline even if the historic development was slightly different. 


2.1. Theoretical Background 

Let t^ = ^ ” ” different points of the Euclidean 

2 

plane H and let z^, z be n real numbers: the observations. As 

1 * n 

in the one-dimensional case, we seek for a smooth function g, called the 
Interpolant , such that 

g(t^) “ i “ X»2 n . (2.1.1) 

In order to define a function Interpolating the data we have first to 

2 

define a domain In ]R containing the data points. There are many possible 
choices and in some cases it will be defined by the problem itself, but while 
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we have the choice we are going to take IR Itself as the domain of defini- 
tion of the interpolant. This will allow us to use the powerful tool of 
Fourier Transform. For other choices, see [ 3 ], [ 29 ] • 

Following the ideas of Theorem 1.1.4, a direct generalization of the 
loss function to be minimized is given by 


J(P) = 


■ 

r ^ 

o 

2 

r 2 ) 

2 

f 9 ) 


3^P 

+ 2 

3 y i 

- 1 - 

iJi 

' 2 

9x^ 


[8x3y^ 


„ 2 
.9y ) 


( 2 . 1 . 2 ) 


In order to set properly the minimization problem we have to introduce 
the following function space: 


where a = 


r^- 2,2 

D L 


= {p: I D%eL^(]R^), |a|=2} 

(aj^,a2>, |a| = 

2 2 

is not equal to H (]R ) . Indeed we have 

22 -222 
H^(]R ) c D L (H ) 


and 


-2 2 
D L 


ts-2. 2 

D L 


is indeed a space of continuous functions (see [ 13 ] , 
is a Hilbert space with the norm 


[J(P)]' 


and 


pep 


(2.1.3) 


(2.1.4) 
[ 14 ]) 


(2.1.5) 


and is the space of polynomials of degree 1 . 

Theorem 2.1.1 . If (t^}. , contains at least one ?, -unisolvent set, 

i” 1 

-2 2 2 

then there exists one and only one element a belonging to D L (!R ) such 
that 


min 

peD“\®(lR^) 

p(t^)*z^ 

i=l,...,n 


J(a) 


J(y) 


( 2 . 1 . 6 ) 
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o is called the Thin Plate Spline Interpolating z^, z^ at 

t^ t" 

D~2^2 

Proof . — = Is a Hilbert space continuously Inbedded In — = — . The set 


= {V€d"V( 1R^) I y(t^)-z^, 1*1,..., n} 


(2.1.7) 


Is clearly nonvoid, and 


I = I + F, 
z z i 


xs Closed In — r since I + P, Is closed In — 

P, z 1 P, 


I + P = n {5 Az )} + P 
^ 1=1 t^ i 1 


and the Injection from — r Into — r — Is continuous. 


Then there exists a unique element (J e 


-2 2 
D L 


.|2 


minimizing jy| over 


The uniqueness of a comes from the fact that any two elements In d 
are different only In a oolynomlal of degree 1 and from the hypothesis we 
have Iq n P^ * {O}. // 


In order to obtain the characterization of the solution we have to use 
convex analysis. Let 

' 0 y e I 


Xj(y) 


otherwise 


Our problem can then be put In the form 


minimize {J(y) + 
yeD’^L^I^) ^ 
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and the solution a satisfies the Inclusion equation 


0 e 3(J + Xj.)(a) . (2.1.8) 


Usi^g now the subdifferential calculus (see [ 19 ]) we get 

n 

0 - - I X, e , (2.1.9) 

in ^ 

2 

where the X, are coefficients, 6 are Dirac measures and A o is the 

2 

Iterated Laplaclan. Moreover, it is clear that the measure A o must be 
orthogonal to Pj^. Then 

n 1 

I X. P(t^) =0 , p e P (2.1.10) 

i=l 

Now we use the fact that (cf. [ 28 ]) 


A^(r^ log r) = 6 

to conclude that 


a(t) = I XJt-t^I^ log It-t^j + q(t) (2.1.11) 

i=l 


where q is a polynomial of degree 1. 

For a complete proof of these results, see [ 13 ] . 


2.2. The Computation of a 

From (2.1.11) - (2.1.10) we obtain the linear system to be solved in 

X^, •••, ■ 


- 17 - 


I A log 

1-1 


Vi * “ 2^2 


n 

I 

1=1 


I 


i 

1=1 ^ ^ 


I A, - 

1=1 ^ 


= 0 


= 0 


0 , 


+ tt3 - Zj , 


j — 1, •••,]! 


( 2 . 2 . 1 ) 


Or, In matrix form 


KA + Ea » 2 

E^A = 0 


(2.2.2) 


where K Is an nxn gymmetrlc matrix and E is nx3 matrix. 


(K)^j = It^-t^l^ log |t^-t^ 


(2.2.3) 


(E)ii 

(E)i, 


= t 


= t. 


(2.2.4) 


(E) 


13 


Many problems arise. The first one is that if K and E are full matrices 

O 

then the solution of this system requires 0(n ) operations. This was not the 
case for one- dimensional splines where the matrix involved in the calculations 
had a band structure and the system can then be solved in 0(n) operations. 


A second major practical problem is that (2.2.2) is ill-conditioned and a 
preconditioning will be necessary to solve it successfully. Finally, from 
(2.1.11) we see that the evaluation of 0 at t requires to use all the 
coefficients. This tends to be unstable due to the alternation of signs of 
the coefficients specially when n is large (depending upon the computer) . 

The ill conditioning of the system can be virtually eliminated using the 
technique we describe now. 

Let Q, R be the Q-R decomposition of E. That is, let Q be an 
n><n orthogonal matrix and R be a 3^3 upper triangular nonsingular matrix, 
such that 


Q*E = 


R 

0 


(2.2.5) 


And let us partition Q in the following way 


1 * = 


L«!J 


( 2 . 2 . 6 ) 


where Q* is 3xn and is (n-3)xn. 

Let 

A = Q*A 


Q*A 

Lqsa 


.A2 


(2.2.7) 


where A^^ e and A2 ^ lT~^. 


The second equation now becomes; 
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OF P00»^ QUr.Ut'V 

E^h - 0 

e'^qa - 0 

[R^IOJA - 0 

r’^a^ - 0 

- 0 . (2.2.8) 

The last equation Is obtained using the nonslngularlty of R. R Is non- 
nonslngular since t^, ..., t*' contains at least a P^-unlsolvent set. 

VJe now use the first equation 

KQA + Ea ® z 

^1 

[KQ. I KQ2I ^ 

^2 

KQ 2 A 2 + Qj^Ra « 2 . (2.2.9) 

Multiplying now both sides by Qi^ we obtain 

Q5KQ2A2 + Q^Q^j^Ra - Q^z 
or 

Q* K Q2A2 = 

since 

QJQl “ 0 . 

Or, finally 

AA2 » z (2.2.10) 
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where 

A - Q*KQ2 
z - Q^z 

a can be easily obtained from (2.2.9) 

Q^Ra • z - KQ 2 A 2 
or 

QJQj^Ra - Q*z - QPQ 2 A 2 

a - R"^[Q*z - QJKQ 2 A 2 I . (2.2.11) 

lor an efficient computer evaluation of the products QJz, Q*z, Q*Kp 2 , 
let us simply observe that 


Q*z 

. J 

Q*z 


Q*K Q - 

' QfKQ^ 

QJKQ 2 ■ 
Q$KQ2 . 


Now the computation of Q*z and Q*KQ is easily done since Q, Q* are 
obtained by a product of Householder transformations that is available for the 
computation of such products when using a package like UNPACK (see [ 12 ]). 

In order to solve (2.2.10) we first prove the following. 

Lemma 2.2.1 . The matrix Q|^ K Q 2 is positive definite. 

Proof. The columns of generate E and the columns of Q 2 generate E"*", 
the orthogonal subspace of E. Then Q 2 represents the orthogonal projection 
onto E'^ written in terms of the basis given by the columns of Q 2 * Thus 


21 - 


y - Q 2 X « E'*’ , 

On the other hand the function 

■ I log It-t^l 

1-1 

with y « £■*■ being a spline function and 

0 < J(0)) - y^Ky 

since from (2.1.9) 

n . 

J(cj) " I V. w(t ) , 

i-1 

and J(a>) cannot be zero since u is not a polynomial of degree 1. // 

Now the system (2.2.10) can be solved using Cholesky's factorization. 

As an Interesting fact we might notice that A does not depend on the 
data itself but only on the data points. So if many data taken at the same 
data points are to be processed, the storage of the Cholesky's factorization of 
A allows us to save a lot of computer time. 

For a complete description of this method and a discussion of its numeri- 
cal properties, see [ 21 ], [ 22 ], [ 43 ]. 

As we have already said, the computation of thin plate splines can be 
very expensive if a large number of points is to be used. In the next section 
we develop a practical method to avoid this problem in some cases. Here we 
must also mention that the actual work towards the development of multidimen- 
sional B-splines could be of great help when dealing with a large number of 
points (cf. [ 7 ]» [ 8 ], [ 26 ]). 

I 


( 2 . 2 . 12 ) 


(2.2.13) 



i 
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2.3. A Practical Method for a Large Number of Polt.ta 

In this section we assume that the number of data points is large and 
well distributed within a rectangle [a,b] [c,d] c ]^. The idea (cf. [ 21 ], 

[ 22 ]) is to divide R, the given rectangle* into several smaller rectangles 
overlapping each other and containing a reasonable number of points. We then 
calculate a thin plate spline in each piece and stick the pieces together using 
a partition of unity. 


More precisely, let 

« - ao < < •• 

• < a < b, a < b* < b, < • • • < b ■ b 

and c * Cq < Cj^ < • • • < 

< c, d < dg 

< dj^ < • • • < dj^ ■ d be partitions of 

[a,b] and [c,d] respectively such that 


‘‘l " 

\ 1- 

0 , . . . , m 


d 1 - 

0 , . . . , k 

1 

1 

(2.3.1) 

®i+l ^ 

^ i - 

0, . . . ,in~l 

*^1+1 ^ 

di i - 

0,...,k . 

These properties Is^ly that 


<*i’V 


0 

(Cj.d^) 


0 


U - [a,b] 

i-1 ^ ^ 


k 

U 

J-0 




[c.dl 
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Ik 


ti» * • • 

OF poOi^ QjAjmr 


and 


D k 

U U Ia,,bjx[c dj - [a,b]x[c.d] 

i-1 j-0 ^ ^ J J 


The situation is illustrated in the following figure 
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We now assume that the number of data points included in 
^ij " ^ "reasonable,” This means that: 1) The 

computer time spent in the computation of a thin plate spline using that number 
of points is reasonable according to the particular problem. 2) The number of 
points is sufficiently large to represent the beh»*^ior of the function in 
tliat region. 3) The niimber of points in the intersection of the regions is 
enough to guarantee a smooth transfer from one region to another. 

Of course all these conditions are extremely subjective and the user has 
to be experienced to know how to apply these conditions to his particular 
problem. 
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Once the subdivision of R in R. . has been deternined, let o be 

Xj Xj 

the thin plate spline interpolating the data on the points belonging to R^^ . 
That is, 


where 


J(o .) ■ min J(y) 

ycD‘*L^(]F?) 

y(t*^)-z. 


(2.3,2) 


k«I 


ij 


■ {k«{l n) i t « R^j"Ia^,b^]x[Cj,d^]} 

Thus is going to be a "good" approximation of the unknown function 

on R In order to build an Interpolant over R, we now "stick" the pieces 
Ij • 

together using a partition of unity on the intersection of the rectangles. 

Let w; IR -► ]R be defined by 


w(x) 


(x-l)^(2x+l) 

0 

1 


X « [o,i; 

X > 1 
X < 0 


(2.3.3) 


Then (u < C (R) and 


a)(x) + (1 - w(x) ) ■ 1 


X e R 


(2.3.4) 


Theorem 2.3.1. Let c.. be thin plate splines interpolating the data 

~ A D 

over the rectangles R^ and Rg respectively. 



*0 


b 


0 
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L. 

or 


Let a.„ be defined by 
Ao 


o^(x.y) 


o^(x»y) 


I o g(x,y) 


(x.y) ® ^ ^ 


p(x)c^(x,y) + (l-p(x))ag(x,y) (x,y) ® ^ 


(x,y) c Rg \ 


where 


p(x) = CO 


X - a. 


'’o-“o 


Then 0^ interpolates the data on u R^ and has continuous first 
partial derivatives . 

Proof. It is clear that interpolates the data on R^ u Rg - R^nRg. 

Let t e R^nR^, then 

o^(t^) = p(x^)o^(tS + (1 - p(x^))Og(A 


p(xS«j^ + (1 - p(x ))Zj^ = 


where t^ = (x^,y^) . 


3a 


V/e only need to prove the continuity of a^g and 


AD 


3x 


along X = a^, 


X = Lq. Both cases being similar we prove it along x * a^. We have 


lim a..(x,y) = lim {p(x)c.(x,y) + (l-p(x))o (x,y)} 


x-*a 


+ AB 


x->-a;! 


a^(ao.y) . 


3a 


AB 


Thus o._ is continuous across the liiie x = a». For - 
AB U oX 


we have 


- 26 - C- - 


3a 


AB 


3x 


i^c^A(x,y) 3o (x,y) 

(x,y) - p’(x)a^(x,y) + p(x) — + (l-p(x)) ^ 


- P’(x) 


30g(x,y) 

3x 


9a 


3x 


and the result follows. // 

Using this method we can stick together all the pieces on each horizontal 
strip and then proceed in the same way to stick the strips together. The final 
result is the following 


G(t) 




fx - a 


u 


i+1 


^l"^l+l 


Oi.Ct) + 


1-03 


X - a 


i+1 


‘’i-^l+l 


t ^ Rij \ ^Rij-l ^ Rij+l Ri+lj ^ Ri-lj ^ 


^i+lj 


t € \ ^Rij +1 ^Rij- 1 ^ 


y-c 


0) 


l+I 




CijCt) + 


1-6J 

•n 

0 

1 

d -c , 


L J J+lJj 


°ij+l^*^^ 


"Rij'^ Rij+l \^Ri+ij^Ri-ij 


(2.3.4) 


0) 


X - a 


i+1 


^i"^i+l 


1-03 








1-03 




l'’3‘"3+ljj 




X - a 


i +1 





03 

fy-c 1 
^ .1+1 


d . — c 

/ J 


L J j+iJ 




1-03 


y- c 


±tl 






t e Rij n n n • 


4 
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It has been proved by Duchon [ 14 ] that this interpolant keeps the 
nice convergence properties of a thin plate spline. On the other hand our 
numerical experience shows that It is a very appropriate method to deal with 
a large number of points. 

Remark . In the above formula we have made the convention 

°0,j “ ^ttri-l.j “ ^1,0 “ ^i,k+l ■ ° * 


3. A General Framework 

The variational properties of natural polynomial splines in one dimension 
lead Atteia [ 3 ] to give a general variational formulation for spline 
functions. This idea is to consider three Hilbert spaces X, Y, Z and two 
continuous linear operators T, A from X onto Y and X onto Z respec- 
tively, then given z ^ Z a general spline of interpolation is the solution to 

min ilTyllY 

yex 

A(y)=z 

All the known cases of splines satisfying variational properties can be 
found using this general approach. For details see [19 ]• Nevertheless, 

in some cases it is somewhat difficult or not natural to use some Hilbert 
space to make the splines fit into this model. For that reason Duchon [ 14 ] 

introduced the semi-Hilbert space model where he works with semi-reproducing 
kernels. In this section we are going to give a slightly different formulation 
using reproducing kernels because of its immediate stochastic interpretation. 
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3.1. Reproducing Kernels and Spline Functions 

Let H be a repfoducing kernel Hilbert space (see [ 2 ], [ 27 ]) 

that is, a Hilbert space where the evaluation functionals are continuous. 

H must be a subset of ]R , the space of functions from X into IR. 

H c . (3.1.1) 

Following Aronzajn [ 2 ] this implies the existence of a function 

K, the reproducing kernel, such that 

i) K(*,x) e H any x e 

ii) K(x,y) = K(y,x) 

iii) f(x) “ {f,K(*,x)) any f e 

where (*,*) denotes the scalar product in 

Let X,, .... X be continuous linear 
1* ’ n 

let z e be an arbitrary vector of real 

Interpolating z^ z^ with respect to 

the unique solution to 

min lly || . (3.1.3) 

yeH 

X^(u)=z^ 

i=l n 

The existence and uniqueness of o the solution to (3.1.3) is a result 
of the projection theorem in a Hilbert space and the fact that {yeH | X^(u)=z^} 
is a closed linear manifold in H. 

The special interest of working with a reproducing kernel Hilbert space 
is based in the following result. 


H 

H. 

functionals from H onto IR. And 

numbers. We define the spline 

the functionals X, , .... X as 
1’ ’ n 
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Theorem 3.1.1 . The unique solution to (3.1.3) is given by 


n 

6(x) = I a h (x) 

i«l ^ ^ 

where 

h^(x) * A^(K(»,x)) 


and a^, satisfy the linear system 


(3.1.4) 


(3.1.5) 


n 

y a. X.(h.) = z , i = l,...,n . (3.1.6) 

j=l J ^ J i 

Proof . Let A: H -*■ r” be the linear transformation 

A(y) = (X^(p), .... X^(y))^ 

It is clear that 6, the solution to (3.1.3), is the orthogonal projec- 
tion of 0 onto the linear closed manifold 

M = {ycH I X^(u)=z^, i=l,...,n} 

exists and is unique if and only if M is non-void. Moreover, we have 
the orthogonality of the projection; 

{o-a,y) - 0 if A(y) = 0 

In other words, 6 e N(A)'*‘ = R(A’). But 


n 

A' (x, , . . . ,x^) - y x.h. 

1 n jii 1 1 


where h^(t) «= X^(K(»,t)) is the representation of guaranteed by the 
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Rltz theorem. Then there exist a,, a such that 

1 n 




n 


i-1 




But d also belongs to M, hence 


n 


X^(6) 


1=1 J J '* 


This completes the proof. // 

As an example of use of this theorem we could try to get the thin plate 
spl ine . 

Of course here the functionals ..., are given by 


X^(y) = y(t ) 


i “ l^...^n . 


(3.1.7) 


12 3 

Let us assume that {t ,t ,t } is a P^-unisolvent set. And let 
{pi, P 2 iP 2 ^ ® basis of satisfying 

3 


I 


k=l 




(3.1.8) 


-2 2 2 

Then it can be shown that (cf. [ 43 ]) D L (]R ) is a Hilbert space 


with norm 


k-1 


IR' 


2 1 

:2 

f 2 1 

2 

f 2 I 

na 

+ 2 


+ 

IJi 


1 

3x9y 

1 J 




(3.1.9) 


and reproducing kernel 
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K(t,s) 


^ It-si^ log |t-s| - I Pj^(t) Is-t’^l^log |s-t^l 

k*l 


87T 


" "Iri ^ log |t-t^ 


k*l 


(3.1.10) 


j 3 

^ Pi<t)p.(s) |t^-t^l^ log |t^-t^l + I p (t) p (s) . 

® l,j-l J k-1 


Now using the theorem, the thin plate splines can be found solving the 
linear system 


I a K(t,t^) = z 

j=l ^ ^ 


j “ l,...,n 


(3.1.11) 


3.2. The Stochastic Approach 

Let 0) be a random process Indexed by X with covariance 

E[ 0 ) CO ] = K(s,t) (3.2.1) 

S w 

and zero mean. 

2 

Let L (co) be the Hilbert space spanned by co^ , x e x and their limits 
In quadratic mean. l^(co) has the Inner product 

l\i,v] - E[pv] . (3.2.2) 

As In the last section, we denote by H the Hilbert space with the ’repro- 
ducing kernel K . 

It Is a well-known result (ef. [ 10 ]) that H and L (co) are 
Isomorphic with the correspondence 

f(x) - E[cc|j^y] (3.2.3) 


f c H 


y e l-^(co) 
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Let A^ be the linear functional of theprecedlng section. It 

is known that A^^, ...» A^ can be defined on co using the isometry. Thus 
A^(o)) is going to be an element of I (oj) defined by the equation 

h^(x) * E[03 ^A^(w)] . (3.2.4) 


We consider now the problem of finding thebest linear estimator for 

Cl) given A. (co) , ..., A (co) . That is, we want to find 
XX n 


CO = E[co I A, ( 0 )),...,A (co)] 
X X j. n 


(3.2.5) 


If we assume that the process is Gaussian, we can find as the least 

2 

squares predictor of co , l.e., as the projection in L (co) of co onto 

X X 

the linear space spanned by A^(co), ..., A^(co) . 

Thus, S is the solution to 

* X 


minimize 

(a , .. .,a )e]R" 
1 n 


i=l 


- “x 


inio) 


and 


(3.2.6) 


n 

cj * y 6t. A, (co) . (3.2.7) 

* i=l ^ ^ 


In order to find 5^, i = l,...,n, we have to solve the system of linear 
equations given by the orthogonality conditions 


E 


f 

n 

I a.A (co) 
i-1 

V 


CO. 


Aj (co) 


0 


j ■ 1, . . .,n , (3.2.8) 


or, equivalently 
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and using the isometry again 




Finally, using the fact that h^ is the representer of we get 


U 

I 5 X (h ) « h (x) 

i-1 ^ J 


(3.2.9) 


If we now assume that z, , ..., z are realizations of the random 

1* * n 


variables X^^Co)), ...» X^(o)) we get 


where 


I 

i=l ^ 


(3.2.10) 


Z ^ (h^(x), ..., h^(x))^ 


in other words 


CO * z^ Z ^ (h. (x), . ... h (x))^ 
x in 


(3.2.11) 


If we now go back to Theorem 3.1.1 we see that 


6{x) “ w 

X 


Theorem 3.2.1 . Let co be a random Gaussian process Indexed by X with zero 
mean and covariance function K(*,*). Then 


9(x) » E[co^ I X^(co)-z^, i«l,...,n] 

For further results see [ 1 , [ ]» I ]» I 1 


(3.2.12) 
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4. Constrained Interpolation 

In the preceding sections we have considered the problem of Interpolation 

of an unknown function given Its values at n are generally distributed points 

2 

of the Euclidean space IR . The hypothesis made on the function Is generally 

-2 2 2 

one of the kind: "f Is smooth." This meaning, most times f e D L (]R ) 

or something similar, but In all cases we have assumed that f belongs to 

some linear space and have used as Interpolant the mlnlmlzer of a semi-norm 

over the linear manifold of Intorpolants belonging to that linear space. As a 

result, the interpolant Is linear as a function of the data and the algorithms 

necessary to compute It require the solution of matrix problems. Even If this 

Is sufficient for many problems, there are others where some nonlinear constraint 

Is naturally posed, as In the case of data coming from a positive function. 

More generally, we could assume that some property of the function is known In 

-2 2 2 

the form f ^ C where C is a convex closed set of D L (]R ) . The problem 
is to find a "good" interpolant satisfying the nonlinear constraint f « C. 

This problem has recently received much attention In one and two dimensions. 
The interested reader should see the literature [ 18 ], [ 32 ] , [35 ], 

[36 ], and the references therein, 

4.1. The Positive Thin Plate Spline 

For these lectures we would only be interested in positive thin plate 

splines (cf. [ 35 ], [ 36 ]). Let z^, ...» z^ be positive real numbers 

2 

and let C be a compact subset of IR . We define the Positive Thin Plate 
Spline as the unique solution to the minimization problem (cf. Th. 4.1.1). 
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minimize J(u) 
y(t) SO, teC 


(A. 1.1) 


where, as In the preceding sections J Is defined by 


J(y) - 


► 

m m 

2 






+ 2 

IJL. 

2 

J. 


le 


1 a* 

9x9y 

• a 

1 

CM 


We first give an existence and uniqueness theorem. 

Theorem A. 1.1 . Let C be compact and assume that {t^, t^, ..., t*'} contains 
at least a P^-unlsolvent set. Then there exists a unique solution of (A. 1.1). 
Proof . Let Pj^ be the set of polynomials of degree 1, then H • ^ 

Is a Hilbert space continuously Imbedded in — r — . The set 


S - {yeD“V(m^) |y(t^)-z^, 1-1 y(t)so, t*C) 

is clearly non-void and 


S - S + Pj^ c H 


.3R' 


is closed since S + P, Is closed ii 

1 P, 


n 


s + p - <n n 5l'^(Io,“))f + p 

[1-1 t tec ' 


r-1, 




.IR' 


and H Is continuously Imbedded In 


Then there exists a uniqj«i element d e H minimizing J(y) over S 
(J(y) Is a norm in H.) The uniqueness of o comes from the fact thaL 


- 36 - 


% 


n {ycD“^L^(»^) I y(t^)-0, 1-1, ...,n) - 10} , 

since {t^, contains a Pj^-unlsolvent set. // 

It Is also possible to give a general characterization jf the solution 
o. We give this characterization In Th. A.1.2 but we omit the proof. The 
Interested reader Is referred to [ 35 ] . 

Theorem 4.1.2 . Let a be the unique solution of (4.1.1) and define 

A - {t e C 1 o(t)-0} . (4.1.2) 

A Is a closed set. Then there exist a positive measure v with support 
(v) c A and n real numbers ...., X^ such that 

n 

^ - I ^ + P (4.1.3) 

1-1 i ^ 2 

where P ^ ^2 elementary solution of the blharmonlc operator 

K2(r) - -^P^(r^logr) . (4.1.4) 

Here * denotes the product of convolution. // 

This last result Is only of theoretical Importance since It has not been 
possible to use It in a numerical algorithm allowing a to be found. In the 
next section we present a numerical procedure allowing the solution to be found. 

For convergence properties of constrained splines see [ 32 ], ( 35 ] . 

4.2. The Dual Algorithm 

Our aim In this section Is to present a special version of the general 
algorithm by Laurent to solve constrained problems by dual iteration (cf. [ 20 ]). 
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j. r .. 

POOR 


We also mention the special numerical methods used to compute the iterate. 

The idea of the algorithm is to construct a sequence of half -spaces 
and solve in each Iteration the constrained problem. 




min J(y) 

y«D"*L*(lR*) 

y(t l"lf . . . »n 

y«Dk 


(4.2.1) 


The sequence {o^} is then shown to satisfy the following properties 
(cf. [ 20 ]) 

i) J(Ok> ^ ^ J(o) 

11) Dk ^ {y«D‘V(»^) ly(t)20, t«c} 
iii) lim 0 , ■ o 


2 

where the convergence can be taken over H (C). 

The starting point of the sequence is the unconstrained thin plate spline. 
That is, the solution of 


J(o ) - min J(y) 

° y«D'^L^(lR^) 

y(t^)«z^, 1-1 n 


2 

Step 0 . Let b^ • K be a solution of 


(4.2.2) 




min o(t) 
t«C 


(4.2.3) 


0RI6?NAL PAC,; 

OF POOR QUALITY 


If ^ 0 ^ " ^0 algorithm stops. Otherwise set k • 0 

and proceed to step 1. 

Step 1 . Let 

I y(bQ)aO} (4. 2. A) 

and let be the unique solution of 

mlnlmlzf, J(w) . (A. 2. 3) 

y«D‘^L^(lK^) 

y<t^)"2j^, i“l,...,n 

y «D^ 

And set b^ as a solution to 

o^(b.) • min o (t) . (A. 2. 6) 

t«C 


Step 2 . Let k k+1. 

If we are lucky enough Cfj^(bj^) ^ 0 and o - Oj^. Stop. 
Otherwise, let 


Dj^ - Dj^ n {y«D"\^(K^) | y(bj^)sO) 


and let 0, . , 
k+1 


be the solution to 


(A. 2. 7) 


minimize J(y) . (A. 2. 8) 

y«D’^L^(lR^) 

y(t^)-z^, l«l,...,n 


Using now the Kuhn-Tuckcr conditions we know that there exist positive 
k 

Vj, j ■ l,...,nj^ such that if we set 


constants 
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we have 


OR/Qir-A:. r 
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^k+1 


y e D 




k k. 


I vj y(b")-o 

l«i '* 


(4.2.9) 


“ ®in J(P) 

y(t^)*z^, 1=1,..., n 


y e D, 


k+1 


(4.2.10) 


Now define b, as a solution to 
k+1 




“in o (t) 
t^C ^ 


(4.2.11) 


and go to step k. 

The preceding calculations are repeated until a convergence criterion is 
satisfied. Generally this criterion is given in the form 

\+l^\+l^ “ (4.2.12) 

where e > 0 is a tolerance givai by the user. 

Remarks 

(1) The computation of is not as hard as it seems to be since 

o^(bQ) = 0 . (4.2.13) 

Otherwise the constraint y e would not be active and the solution 
might give <^^^^0^ ^ ^ which is Impossible since in that case Oq = 

Thus is Indeed the solution to the simple interpolation problem 


(4.2.14) 


minimize J(u) 

yeD'^L^Cm^) 

y(t^)= 2 ^, 1-1,... ,n 

y(bQ)»0 

And Oj^ is of the form 

o^(t) » I X^|t^-t|^ log It^-tl + alt-bpl^ log |t-bQ| + q(t) 

1=1 

with q « P^. 

Then is obtained solving a linear system. 

Similar considerations apply for the computation of step k. 

(2) From the last remark it appears that each iteration costs the solution 
of one or two linear systems where the matrices which depend on the knots are 
different from one step to the other. This could be very expensive to do unless 
we know something more about the systems to solve. 

Here the key remark is that the difference between the original matrix 
(to compute and the matrix used in the computation of is only given 

by a rank 1 perturbation. Then we can use this fact to solve the system in 
step k using the solution in step 0 and the Cholesky factors obtained there. 
(For details see [ 35 ]) • 

In a recent Ph.D. thesis presented at the University of Wisconsin, M. Villa- 
lobos has used a different approach to give an approximate solution to (4.1.1). 
His approach basically consists in the discretization of the constraints, imposing 
a regular grid over C and replacing ij(t) ^0, t e c by 

Vigh SO, g'" ^ G , 

G: grid over C. (For details, see ( 36 ]). 
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5. Thin Plate Smoothing Splines 

In the preceding sections we have considered the solution of some Inter- 
polation problems In two dimensions. In all the cases we have assumed that 

the data are exact and we wanted to obtain a smooth surface passing through the 
1 3 

points (t ,z^) e ]R . Our aim in this section is to present the case of 
noisy data. Thus, in the first part we present the theoretical results giving 
a characterization of the thin plate smoothing spline and in the second part 
we give the numerical methods used in the computation of such splines. 


5.1. Existence, Uniqueness and Characterization 
As in Section 1.2, we assume the model 

z^ = f(t^) + i “ 1,2, ...,n (5.1.1) 

-2 2 2 . 

where f e d L (R ; is unknown and the ^re supposed to be realiza- 

tions of independent identically distributed Gaussian variables. 

In order to approximate f, we introduce the thin plate smoothing spline 
of parameter X > 0 as the unique solution to (cf. [13 ], [14 ]) 


minimize 



(5.1.2) 


As in the one-dimensional case, X controls the tradeoff between the 
smoothness of the solution measured by J(y) and the approximation to the 


1 

n 


1=1 


(y(t^) - z^)^ 


k 


data measured by 


L’vy 


- 42 - 


OF POOR Qui, ;.!,Y 

In order to prove the existence and unlqueuess of o , , the solution 

n, A 

to (5.1.2), let us first define the nxn symmetric semi-positive definite matrix 
as 

y “ min J(y) . (5.1.3) 

yeD’2L®(m^) 

u(t^)=yj^, i-l,...,n 

y^fiy represents the value of J a; the thin plate spline inter- 

polating y^ at t^, i = l,...,n. Given that y ■* is linear, the trans- 
formation y is quadratic hence th<; existence of 

Now we write (5.1.2) in the fol] owing way 


min min •|xj(y) + - T (y(t^)-z.)' 

yelp yeD’^L^(m^) [ "i=l ^ 

y(t^)»y. , i=l n 


= min 




min UJ(y) + 7 j (y^-zJ‘ 
yeD- 2 L=^(lR 2 ) [ "i=l ^ ^ 

[y(t^)=y^, i=l n 


min -s ■< 
ye]^* 


X min 


yeD"^L^(]R^) 


J(y) 


1 " 

n ^ ^^1 

i*l 




[ y(t )=y^,i=l nj 


min 

ye^ 


jxy'Ky + 


(5.1.4) 


The solution y to this problem is easily found by differentiating and 
setting the gradient equal to zero . We obtain 


(nXn + I)y * z 


(5.1.5) 
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And o , is obtained interpolating y at t^, i - 1 n using 

11 f A X 

a thin plate spline. As in the one-dimensional case, (5.1.5) is only of theore- 
tical use as we shall see later. A more practical result is given by; 

Theorem 5.1.1 . Let {t^, ...» t"} contain at least one Pj^-unlsolvent set. 

Then there exists one and only one solution , to (5.1.2). Moreover, there 

exist n+3 real numbers X, , ..., X ; a, , a_, a, such that o . has the 

1 n 1 2 3 n,X 

expression 

^ 12 i 

^(t) = I X^|t-t I log It-t I + o^tj^ + a^t2 + <^3 ( 5 . 1 . 6 ) 


where t = (tj^,t^) • 

The coefficients 

n 

STmXX. + ^ 

j=l 

n 

I 

3-1 

n 

I 

j-i 

n 

I 



X. K(t^-t^) + u,tj + 

3 11 

X.t^ - 0 


X tJ 
j 2 


0 


= 0 


are the solutions of the system 

^ 2^2 ^ ^ " 


(5.1.7) 


Proof . See [ 13 ], [ 14 ]. 

In the next section we use this formulation to compute the smoothing spline. 

Now we are interested in obtaining an interpretation of a . as a filtering 

n ^ A 

process. To do that, let ~ ^2 “ * * * “ ^n eigenvalues of in 

increasing order and let Q be the orthogonal matrix diagonalizing ^2. Then 
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OF POOR QUhLIh f 


where 


n - Q ^ Q* 


Z “ dlag (Pj^, ..., p^) 


Then (5.1.5) becomes 


( 5 . 1 . 8 ) 


where 


Hence 


y * 


(I + nXZ)‘^£ 


Q*y 

Q*z 


( 5 . 1 . 9 ) 


1+nXp^ ‘i 


( 5 . 1 . 10 ) 


But it has been observed, and proved in some cases [ 30 ), [ 39 ] that 


then 


Pi "IT 


^i r?ai^"i 


( 5 . 1 . 11 ) 


which turns to be again a low-pass Butterworth type filter. For details see 
I 35 ], [ 39 ], [ 43 1. 


5.2. The Computation of A and ^ 

As in the one- dimensional case, the estimation of X from the data is a 
very delicate problem. One of the methods that is becoming the most popular 
to solve this problem is the Generalized Cross-Validation which amounts to 
minimizing the GCV function 




If poon Q;j/'.LiT^ 


y/vx<^^^ 

fl - itr.(A(X))1^ 


(5.2.1) 


where A{X) is given by 


A(X) - (I + nXn)' 


(5.2.2) 


Thus, we need an algorithm to compute 0 . and one to compute 

n ^ A 

tr. (A(X)). We begin with the computation of x' 

As we have seen In Th. 5.1.1, the computation of 0^- involves the solu- 
tion of a linear system in the n+3 coefficients X^^, ..., X^; a^t cx^. 

The linear system to be solved is given by equations (5.1.7) or, using the 
notation of Section 2.2, in matrix form 


(SiTnX + K)A + Ea * 2 


(5.2.3) 


Thus, the linear system to be solved has the same form as (2.2.2) and 
using the same technique, let Q, R be the QR decomposition of E 


where Q* is nxn and can be decomposed as 


and Q* is 3xn, Q* is (n-3)xn. 


w 


1 
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Now setting 





‘ QJA ‘ 

A « 

1 

1 CM 




Q*A 


the equations become 


0 


Q*(8TTnXI + K)Q2A2 • Q^z 


Ra = Q*z-Q*KQ2A2 


(5.2.4) 


This linear system can be easily solved using Cholesky's decomposition 
and all the remarks given In the case of Interpolation apply to this case. 

The computation of tr. (A(X)) requires additional work. In [ 35 ] 
we proved that In some cases the eigenvalues of 0 . can be well approximated 
by those of a fourth-order differential operator, but In two dimensions this 
is too complicated to calculate and a direct evaluation of the eigenvalues 
Is necessary. We first obtain an expression for to do this we consider 

the thin plate spline Interpolating y •* (y^,...,y^). Thus its coefficients 
'^1* ^n’ ^1* ^2’ S satisfy the system 

KA + E6 = y 

e’’’a = 0 


and using the notation 


<j 

■K.H 

O' 

1 


— 1 



. ^2 J 


A 


Q*A 
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we get 




(5.2.5) 


0 


Now we observe that from Th. 2.1.1 we have 


J(o ) - AKA 

y 


a"^q* k qa 


(5.2.6) 


A2Q5KQ2A2 


Using now (5.2.5) we finally get 


J(Oy) 


and 


y^Q2(Q*KQ2)'^Q^ 


Q2(Q*KQ2)'^Q* 


(5.2.7) 


Finally, we observe that the eigenvalues of fi are those of Q^Kp 2 
except for the first three that are equal to zero, then 


tr. (I+nXJ^)"^ * 3 + tr. (I + nXq*KQ 2 )"^ 


(5.2.8) 


Thus, if we want to compute tr. (A(X)) many times, as is necessary when 
minimizing the GCV function, we first compute the eigenvalues of Q^KQ 2 and 
then use them to compute tr. (A(X)). For more details and numerical experiences 
the reader is referred to [ 35 ], I 36 [ ^3 ] . 
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6. Non-varlational Techniques 

As we have pointed out before, the methods we have exposed here are 
essentially generalizations of splines In one dimension using the variational 
properties of those functions. We have chosen this way because of Its close 
relationship with stochastic estimation and because It appears to be at 
present the only technique ised for scattered d^ta (when referring to splines). 
Many other generalizations are possible, and we will briefly describe two of them 
in these sections. For a detailed discussion see the references below. 


6.1. Tensor Product Splines 

The first generalization of one- dimensional splines to many dimensions has 


been given using tensor products. More precisely, consider a square [a,b]x[c,d] 
2 


c ■ c, < c- < 


in 3R and two partitions, a ■ a^ < 32 < < a^ ■ b; v. - > ^2 

of [a,b] and [c,d] respectively. Also let 1 ■ 1, ..., n+1 and 

2 

, j “ 0, ..., nH-1; the B-spllnes based on the knots 
{a^,a^,a^,a^,a 2 ,a^,. ..,a^,aj^,a^,a^} and {c^^ , c^ , Cj^, c^^, C 2 , c^, . . . , c^, c^,Cjj^,Cj^} 
respectively . 

^ n+1 

As we have already said (cf. Section 1.1), ® basis for the 

space of twice differentiable piecewise cubic polynomials based on the knots 

( a, , a)* 

■i- n 

We now define a space of tensor product splines S on the grid (a^,Cj), 

1 * l,...,n; j B l,...,m as the linear space spanned by the tensor product 
, lin+l f lim+l 

of the basis j"0 ' s e s If and only if 

n+1 m+1 


< c * d 


m 


1-1 i-O ^ ^ 


i 


( 6 . 1 . 1 ) 
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for (x,y) « [atb]^[c,d] and « ]R for 1 “ 0, n+1; J ■ 0, mfl. 

Or, In other words. 


n+1 ro+1 

I I “ 

1-0 j-0 


Ij i J 


( 6 . 1 . 2 ) 


It is well known (cf . t 1 ], [ 9 ], [ 26 1) that S is a linear 

space of dimension (n+2) ^ (m+2) . Thus, the complete determination of s from 
Interpolation conditions requires to give 2(n+mf2) additional boundary condi- 
tions. A usual choice for these conditions is: 

Hermite type 


3s 

3x 

(aj.Cj) 

- 

given constant , 

j - 

1 > • 

• * y 2H 

3s 

3x 


m 

given constant , 

J - 

1,. 

• • y ID 

3s 

3y 


- 

given constant , 

i - 

1 , • 

* • y Tl 

3y 


m 

given constant , 

i - 

1.. 

■ • y XI 

3^s 

3x3y 


- 

constant 




d 8 

3x3y 


- 

constant 




3^ 

3x3, 


- 

constant 




3 s 
3x3y 

(a ,c ) 
n’ m 

m 

constant 





As in the one-dimensional case these conditions produce optimal convergence 


rates but as it was already pointed out, they require some additional informs- 


tlon on the function that is not always available. (For other choices see 
[ 9 1 . [ 26 ])• 


The functions in the space S are piecewise bicubic. This means 
that in each subre ^tangle function Is a bicubic 

polynomial, l.e., a function of the form 


2 3 

p(x,y) - qQ(x) + yqj^(x) + y q2(x) + y q3(x) 


(6,1.4) 


where q^, q2» q^ cubic polynomials. 


Moreover, this function s ^ S has continuous derivatives 


5s 

dx 


9 s 


i!s 


9 s 


is 

9y ' 


ay 


2 • 


Finally, It Is Interesting to say that the solution of the linear system 
in a^j, j - 0, ...,mfl, 1 » 0, ...,n+l, can be efficiently performed using 

very specialized techniques (cf. [ 1 ], [ 9 ), [ 26 D • For further 

details the reader Is referred to the extensive literature In this area. 


6.2. Multidimensional B-spllnes 

A ctore recent approach to the problem of multidimensional data has been 
given by DeVore, Dahmen, Micchelli and others [7 ], [8 ], [ 26 ] . 

This new approach Is attracting much attention recently, mainly from potential 
users in finite elements codes. 

The idea of this method Is to generalize the B-spllnes to several dimen- 
sions. The idea Is Interesting because the two main problems arising with thin 
plate splines, that is, full system of equations and instability In the evaluation, 
could be solved when using a local support basis. Unfortxinately, until now 
the theory of multivariate B-spllnes seems too complicated to be used in 
practice. For details the reader Is referred to [ 7 ], [ 8 ] and the 
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references therein. Here we only give the basic idea of the construction. 

To do this we go back to the simplest one- dimensional case, the P.-splines of 
degree 1. *^nese are defined as 


X2» x^; t) 


(t-x) 


( 6 . 2 . 1 ) 


where we have made explicit the dependence of the B-spllne on x^^, 
the knots of the spline. A tjpical plot of M is given below 



Of course the use of (6.2.1) to define the multi-dimensional B-spllne 

would require the definition of a multivariate divided difference. Thus, this 

way seems difficult. However we can still generalize (6.2.1) if we observe 

that we can give a geometric interpretation of this definition. To do that, 

2 

let (Xj^.y^^), (x 2 ,y 2 )t (x^.y^) be three points in general position in IR 

(i.e., they form a triangle) 

<4 4 


2 2 
Now for given x ^ JR define fl(x) as the volume In IR of the set 


! (Jt.y) « T) 


( 6 . 2 . 2 ) 


where 


Then 


T - Triangle with vertices (x^.y^), (X 2 ,y 2 )t (x 2 ,y 3 > 
M(x) - voljj^i(Ij^) . (6.2.3) 


It is clear that M has the same form of M(x^,X 2 ,X 3 ; *) except maybe 
for its maximum value. Thus (6.2.3) is an alternate definition for M. Clearly 
this function does not depend on (yj^.y 2 »y 3 ) only on the volume of T. 

S 

The generalization to IR Is now clear. Let 0 be a unit volume simplex 

S+k S 

in IR and define for x « H the function 

M^(x) - I (x,u)«o} . (6. 2. A) 

This function is a "smooth" piecewise polynomial with support given by the 

S 

convex hull of the projection of o in ® . 

Most of the classical results on one-dimensional B-spllnes can be obtained 
with these new functions. For example: recurrence relations, integral forms, 
etc. See I 7 ]. [ 8 ], [ 26 ]. 


6.3. Conclusion and Comments 

Many other subjects might have been Included in a coiqilete treatment of 
two-dimensional or multidimensional splines. The subject is an active area of 
research and the references given below are only a sample of the extensive 
literature in the area. The Interested reader should consult the symposium 
on multivariate approximation and related research meetings to have a better 


idea of the most recent results. 
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